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Abstract — Underground  imaging  of  dielectric  and  conductive 
anomalies  performed  using  ground  penetrating  radars  (GPRs) 
requires  expensive  wideband  systems  to  increase  the  resolution. 
The  advent  of  tomographic  principles  in  multi-monostatic  GPRs 
dramatically  improved  the  imaging  capabilities  and  suggested 
the  possibility  of  reducing  the  bandwidth  of  the  probing 
waveform.  In  this  work  we  propose  to  extend  the  tomographic 
principles  to  the  case  of  below-ground  distributed  sensing,  thus 
taking  advantage  of  the  geometric  diversity.  We  show  that,  by 
using  geometric  diversity,  the  frequency  content  required  to 
image  below-ground  targets  is  drastically  reduced  to  virtually  a 
single  monochromatic  signal,  thus  achieving  full  spectral 
dominance  in  the  waveform  design. 

I.  Introduction 

Presently,  a  prevalent  approach  to  detect,  locate,  trace  and 
image  underground  objects  is  by  using  Ground  Penetrating 
Radar  [  1  ]-[5].  Depth  of  penetration  is  augmented  by  logging 
data  through  boreholes.  Image  resolution  is  ameliorated  by 
using  tomographic  principles  applied  to  the  received  data  [6]- 
[8].  However,  all  GPR  systems  need  to  use  short  pulses  (i.e. 
high  frequency  bandwidth)  to  increase  the  information 
concerning  the  targets  via  frequency  diversity. 

The  use  of  high  bandwidth  leads  to  several  issues.  First, 
the  signal  to  noise  ratio  (SNR)  decreases  with  an  increase  in 
the  spectral  content  of  the  probing  wavefield.  Second,  the 
electromagnetic  spectrum  available  for  military  and  civil 
applications  is  continuously  being  eroded  due  to  the 
tremendous  demand  of  wireless  applications.  Furthermore, 
unintentional  (e.g.  broadcasting  stations)  or  intentional  (e.g. 
jammers)  man-made  interferences  can  reduce  the  available 
spectrum.  If  larger  bandwidth  is  required,  EMI,  EMC  and 
intermodulation  effects  become  difficult  tasks  to  be  tackled 
and  solved.  Third,  a  wideband  system  can  be  extremely 
bulky,  delicate  and  expensive.  This  problem  is  accentuated 
when  the  system  is  designed  to  work  at  lower  frequencies:  it 
is  impractical  to  generate  wrell-designed  short  pulses  in  the 
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HF  frequency  range,  where  antennas  are  fundamentally 
electrically  small.  Additionally,  the  shorter  the  pulse  is,  the 
greater  is  the  difficulty  and  cost  in  properly  sampling  the 
received  signal.  Fourth,  the  conductivity  and  dielectric 
permittivity  of  the  ground  varies  with  the  frequency:  if  a 
wideband  signal  is  sent  into  the  ground,  frequency  dispersion 
is  likely  to  occur,  thus  broadening  the  pulse  support  and 
reducing  the  achievable  resolution.  Additional  non-bandwidth 
related  factors  limit  the  efficiency  of  common  GPRs.  For 
example,  the  resolution  in  azimuth  depends  on  the 
beamwidth:  at  HF  frequencies  it  is  impractical  to  create 
pencil  beams,  therefore  azimuth  resolution  must  be  reduced 
by  using  other  techniques.  Moreover,  common  GPR  suffer 
from  the  blind  region  effect  in  which  the  receiver  is  idle  until 
the  transmitter  completes  the  transmission  of  the  pulse.  This 
problem  can  be  solved  by  invoking  suitable  modulation 
techniques,  at  the  penalty  of  increased  complexity  and  cost  of 
the  system.  Finally,  when  a  target  is  not  parallel  to  the 
surface,  Euler’s  law  suggests  that  reflected  energy  is  not 
mainly  back-propagating  to  the  GPR  receiver,  and  thus  its 
detection  may  be  compromised. 

We  propose  a  methodology  [9]-[  14],  named  RF 
tomography,  that  addresses  these  open  issues,  and  we 
demonstrate  how  RF  tomography  may  outperform  current 
state-of-the-art  GPR  technology.  RF  tomography  requires  a 
set  of  low-cost  penetrable  transponders  arbitrarily  deployed 
above  the  ground  (see  Fig.  1).  Transmitters  send  a  waveform 
into  the  ground,  scatterers  re-irradiate  power  toward 
receivers,  which  log  data  and  relay  the  retrieved  information 
to  a  base  station.  The  novelty  of  this  approach  is  the  use  of 
multiple  transmitters  (i.e.  view  diversity)  and  multiple 
receivers  (i.e.  observation  diversity),  besides  frequency, 
polarization  and  antenna  pattern  diversity.  The  accrual  of 
geometric  diversity  facilitates  the  waveform  frequency 
content  reduction,  and  in  the  limiting  case  a  set  of  discrete 
monochromatic  signals  yield  actionable  below-ground 
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reconstructed  images.  In  this  work  we  focus  on  the  inversion 
algorithms  necessary  to  process  the  signals  measured  by  the 
multitude  of  ground  sensors. 


Figure  1 :  RF  Tomography.  Transmitters  send  power  into  the  ground. 
Receivers  collect  the  scattered  field  and  send  this  information  to  the  main 
station. 


II.  Forward  Model 

We  describe  the  forward  model  of  an  RF  tomographic 
system  by  considering  the  3D  geometry  depicted  in  Fig.  2. 
The  host  medium  (i.e.  the  earth)  is  modeled  as  an 
homogeneous  medium  with  relative  dielectric  permittivity 
£D  ,  conductivity  (JD ,  and  magnetic  permeability  //() .  The 
targets  are  assumed  to  reside  in  the  investigation  domain  D. 
The  sources  are  N  electrically  small  dipoles  (of  length  A/' ) 
or  loops  (of  area  A1 )  fed  with  current  /' ,  and  located  at 
position  r'n  (view  diversity).  For  each  transmitting  antenna, 
the  scattered  field  Es  is  collected  by  M  receivers 
(observation  diversity),  located  at  points  in  space.  For 
simplicity,  a  single  operating  frequency  /  is  adopted. 

We  assume  the  relative  dielectric  permittivity  profile 
£t .  (r')  and  the  conductivity  profile  <T(r')  inside  the 

investigation  domain  D  as  unknowns  of  the  problem. 

Accordingly,  the  inverse  problem  is  recast  in  terms  of  the 
unknown  permittivity  contrast  function: 


^(r')  =  £,(r')-£0  +  7 


g(  r')-go 
27tf£u 


(1) 


In  this  way,  the  wave  number  inside  D  can  be  expressed  as: 


k-  (r')  =  or n0£0£r  (r  ')  + jcofi0cr(r ') 
=  k2D  +  k2a£s(  r') 


(2) 


kD  =  6)J{10£0£d  +  /  O) 

k0=O)JJ^£u 


(3) 


The  function  in  (1)  accounts  for  the  difference  between  the 
unknown  dielectric  permittivity  of  the  object  and  that  of  the 
host  medium. 

For  each  point  r'  in  region  D ,  the  vector  wave  equation 
holds: 


VxVxE(r')=[^+fe(r')]E(r').  (4) 

The  scattered  wave  in  a  point  rg  D  that  is  solution  of  (4) 
can  be  written  in  terms  of  integral  equation  of  the  dyadic 
Green's  function: 

Ei(r)  =  *.2JfJi(r,r’)-E(rV,(r')**.  (5) 

D 

where  E(r ')  is  the  total  field  in  the  investigation  domain  D, 

given  as  the  superposition  of  the  incident  field  E1  (r‘)  (i.e. 
the  field  in  the  investigated  area  when  objects  are  absent)  and 
the  field  E's  (r) ,  scattered  by  the  targets. 

As  it  is  well  known,  the  inverse  scattering  problem  in  (5)  is 
non-linear.  Nevertheless,  it  can  be  recast  to  a  linear  problem 
by  means  of  the  Bom  approximation  (BA).  Under  BA,  the 
total  field  inside  the  integrand  of  (5)  can  be  approximated  by 
the  known  incident  field  [  1 5]-[  1 8],  i.e.: 

E 9  ( r )  =  A:02  JJJ  G  ( r , r  ■)  -  E'  ( r  ■)  ^  ( r ')  • .  (6) 

D 

Therefore,  the  inverse  problem  at  hand  is  cast  as  the  inversion 
of  the  linear  integral  equation  connecting  the  permittivity 
contrast  function  to  the  scattered  field  data. 

The  use  of  BA  can  be  justified  by  considering  that: 

•  The  targets  of  interest  are  isolated,  limited  in  number  and 
embedded  in  a  lossy  medium.  Therefore,  mutual  interaction 
(a  phenomenon  ignored  by  BA)  between  targets  can  be 
assumed  negligible. 

•  In  general,  the  inhomogeneities  of  the  soil  are  electrically 
small,  and  their  conductivity  remains  low.  Therefore,  their 
scattered  fields  are  insignificant  compared  with  the  RF 
signal  re-irradiated  by  our  targets  of  interest. 

•  Our  goal  is  to  detect,  localize  and  approximately  determine 
the  geometry  of  the  targets.  Toward  this  objective  it  has 
been  shown  how  BA  based  inversion  algorithms  are  able  to 
work  with  strong  scattering  objects,  provided  that  no 
quantitative  description  of  the  dielectric  permittivity  in  D  is 
required. 
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The  Green's  function  depends  upon  the  geometry  that  is 
considered  (half-space,  two-dimensional,  full  space).  In  this 
work,  we  utilize  the  three  dimensional  dyadic  Green’s 
function  G  for  a  homogeneous  medium  with  the  same 
properties  of  the  earth.  This  assumption  is  reasonable  because 
the  sensors  are  deployed  at  the  air/ground  interface,  and  the 
frequencies  involved  are  relatively  low.  Accordingly,  we 
have  [19]: 


G(r,r’)  = 


1  + 


VV 


=  ki 


e'  ° 


fcnlr-rl 


4k  r-rl 


(7) 


The  operator  VV  in  (7),  which  is  responsible  for 
depolarization  and  is  useful  for  near  field  sensing,  can  be 

generally  neglected  for  £0|r  —  r'|»l,  i.e.  when  sensors 

are  located  in  the  far  zone  with  respect  to  targets. 

The  incident  field,  i.e.  the  field  radiated  in  the 
homogeneous  medium  from  a  point  source  located  at  position 

rln  ,  is  given  by: 


E:,(r',r,;)  =  OG(r',r:)-a:„ 


(8) 


III.  Inversion  Procedures 

A.  Tikhonov  Regularization 

A  way  to  compute  L-1  is  to  perform  a  numerical  inversion 
of  L  [20].  Let  us  collect  the  sampled  field  data  in  an  ordered 
NM  vector  ES  =  \eS  (  v'n ,  r'  )  j  ,  and  discretize  the  domain 


where  Q  =  jO/i^Al1 11  for  an  electrically  small  dipole,  or 

Q  =  —  jO)/J()A' I'  for  an  electrically  small  loop,  is  the 
(electric  or  magnetic)  dipole  moment  direction. 

Additionally,  the  field  received  by  a  dipole  or  loop  with 
moment  direction  a^f  positioned  at  due  to  an  equivalent 

(in  terms  of  E* )  current  distribution  defined  inside  the 
investigation  domain  D  can  be  expressed  as  [28]: 


region  D  in  K  voxels,  each  one  located  at  position  rA. ' :  the 
contrast  dielectric  permittivity  can  be  embodied  in  a  column 
vector  £$  =  (rA. ')}  of  length  K ,  and  it  represents  the  set 

of  unknown  parameters.  After  this  discretization,  eq.  (10)  can 
be  rewritten  in  a  matrix  form: 

Es  =  L£>-.  (ID 


£s(r>:Mox 

JJJ - G  (r;; , r ')  -  E'  ( r r' ) ^  ( r >/r ’  • 


(9) 


Substituting  (8)  in  (9)  we  obtain  the  scalar  forward  model  of 
RF  tomography: 


Es  (r>;)  =  L(^(r’))  =  0Ao2x 


From  a  mathematical  point  of  view,  the  problem  of  finding 
the  contrast  function  is  to  perform  the  inverse  of  the  linear 
operator  L  connecting  the  unknown  contrast  function  and  the 
scattered  field  data. 


where  L  now  is  a  matrix  with  dimensions  NM  X  K  . 

The  problem  is  then  to  invert  the  relation  (11).  Due  to  the 
independent  set  of  measurements,  L  is  theoretically  full  rank, 
but  is  often  severely  ill-conditioned.  This  leads  to  severe 
artifacts  in  the  reconstruction  process,  particularly 
exacerbated  when  noise  (thermal,  external,  quantization)  or 
clutter  is  impinging  the  receivers. 

A  common  way  to  quantify  the  behavior  of  L  is  by 
inspection  of  its  condition  number  K .  For  the  operator  L  it 
is  quite  common  to  obtain  typical  values  of  K  above  106 . 

An  efficient  method  to  perform  an  inverse  of  a  very  ill 
conditioned  matrix  is  by  using  the  Tikhonov  regularization 
procedure.  In  this  way,  the  contrast  dielectric  permittivity  can 
be  estimated: 


£  =  (L"L  +  ^I)"'l "Es  (12) 
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Where  L.'1  denotes  the  adjoint  of  L ,  and  j5  is  the 
regularization  parameter  in  the  Tikhonov  sense,  that  needs  to 
be  appropriately  selected.  The  advantage  of  this  approach  is 
its  remarkable  performance  in  generating  meaningful  images, 
even  when  the  number  of  sensors  is  limited.  Unfortunately,  a 
proper  choice  of  ft  may  be  a  difficult  task,  and  often  it  is 
necessary  to  seek  for  a  constrained  optimization  solution  of 
P  before  a  meaningful,  sharp  and  low  blurred  image  is 
reconstructed.  This  implies  a  (computationally  expensive) 
matrix  inversion  for  each  attempt  may  be  necessary. 

B.  Fourier  Approach 

In  practical  scenarios,  where  real-time  processing  is  critical, 
or  when  the  relay  to  a  base-station  is  impeded,  it  is  pivotal  to 
derive  an  inversion  strategy  that  privileges  speed  vs. 
accuracy.  This  priority  is  emphasized  by  the  current  system 
technology,  which  is  widely  implementing  FFT  routines  to 
accelerate  image  reconstructions.  Therefore,  we  propose  an 
approach  that  takes  advantage  of  the  Fourier  relation  arising 
between  scattered  field  and  object  shape,  as  discussed  in 
literature  under  the  topic  of  diffraction  tomography. 

In  fact,  if  targets  and  sensors  are  distant  enough  so  that  the 
propagating  wave  is  TEM  (normally  occurring  when  the 
fields  are  primarily  propagating  as  1  /  r ),  then  the  forward 
model  can  be  expressed  as  follows  below. 

We  define  the  unit  norm  direction  of  propagation  vectors  as: 

1,',  =  -x  sin  ff„  cos  <p'n  -  y  sin  6[  sin  ^  -  z  cos  (13) 
K,  =  *  sin  O',,,  cos  <p'm  +  y  sin  0rm  sin  <p'm  +  z  cos  ffm .  (14) 


Using  the  paraxial  approximation,  the  transmitting  Green’s 
function  at  the  generic  position  r'  inside  region  D  can  be 
simplified  as: 


,  ,  ,  exp [+jkDr'a )  exp(+y£0i;  •  r ') 

(r',r  )  =  - - L,  (15) 


4  7Tr[ 


while  the  receiving  Green's  function  can  be  expressed  as: 


4  jtrr 


Therefore,  for  a  pair  of  transmitters  and  receivers,  the 
scatteres  field  can  be  rewritten  as: 


E*( r'  r ')=&'-<  qSM'K 

[  \67T2~'~'  L 


X 


r  r 

>i  in 


•(17) 


JJJ  %  (r ')  exp  (+jkD  [  1'  -  Ym  ]-r')«/r' 

D 

Eq.  (17)  is  a  simplified  and  continuous  version  of  (10). 

The  quantity  kD[\'n—  can  be  represented  by  a  3D 
vector: 


k . =M 


(18) 


Eq.  (17)  can  be  rewritten  as  [33]: 


E\K^=^QeMC+K]X 

K  mn>  16  TTrlr^ 


IfI^(r')exP(+/k™»  r'l'/r' 


(19) 


It  is  useful  to  consider  a  normalized  version  of  (19): 


F(k  )-  [6*2«' 

E  [K"’~Qkl a'-a'  e 


EUk 


=  JJJ (r  ) exP (+yk„„  • r ') dr ' 


(20) 


This  result  can  be  interpreted  in  the  following  way:  each 
collected  sample  ES  (kj)m )  returns  the  value  of  the  kw;f 
spectral  component  of  the  contrast  function 
Theoretically,  if  we  have  enough  samples  to  fully  populate 
the  spectral  representation  of  (r'),  the  discrete  function 

ES  (km/J)  in  the  limit  can  be  approximated  as  a  continuous 
function  Es  (k) ,  and  (20)  can  be  interpreted  as  a  3D  inverse 

Fourier  transform  of  the  permittivity  contrast  function. 
Therefore,  we  can  reconstruct  an  image  of  the  underground 
by  direct  Fourier  transform  eq.  (19),  i.e.: 

^(r')-}jj^S(k)exP(--/kr')</k-  (2D 

K 


where  the  domain  of  integration  K  is  the  support  of 
Es  (k)  .  By  inspection  of  (18),  we  conclude  that  when  the 


sensors  completely  encircle  the  target,  K  is  a  sphere  of 
radius  2 kD,  meaning  that  the  available  information  of  the 
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spectral  content  of  £^(r')  is  limited  up  to  the  spectral 
component  2 kD .  Therefore,  the  reconstructed  image  of  the 

contrast  function  will  be  a  low-pass  filtered  version  of  the 
true  image. 

In  the  real  scenario  where  a  finite  number  of  sensors  are 
deployed,  three  factors  affect  the  resolution  (leading  to 
blurring  and  artifacts):  the  invalidity  of  paraxial 
approximation,  the  non-uniform  sampling,  the  sparse  data  set, 
and  the  attenuation  constant. 

In  this  work,  we  consider  the  attenuation  to  be  negligible, 
so  that  the  wave  number  in  (21)  remains  a  real  quantity,  and 
FFT  can  be  applied. 

Paraxial  approximation  holds  when  the  angle 
0max  between  the  ray  passing  through  the  origin  and  the  ray 

intersecting  the  boundary  of  the  region  D  is  negligible.  This 
angle  can  be  computed  using: 


KyAa)  = 


>-T J— 

AA'„„ 

0  otherwise 


The  total  weighting  factor  is: 


A'  A/ 


w . 


mi  u 
mn 


n=\  #n= 1 


The  estimated  value  in  (w,  v,  w)  is  therefore: 


E  ( 1 « » ■ v> 1 w) = Z  Z  w’’»'E  ( k ' 


n=l  m=\ 


(25) 


(26) 


(27) 


0_ . =  max 


max 

t'gD 


tan 


-i 


,  (22) 


where  /  represents  any  transmitter  or  receiver.  Blur  reduction 
is  accomplished  by  segmenting  the  region  D  into  smaller 
analysis  regions  (where  0ttaK  remains  small  within  the  sub- 
region)  and  by  considering  an  inverse  problem  (i.e.  smaller 
FFT)  for  each  sub-region.  Then,  the  resulting  sub-images  are 
concatenated  to  form  the  final  image. 

The  non-uniform  to  uniform  grid  transformation  can  be 
accomplished  using  Tri-Linear  interpolation.  Let  us  define  a 
uniform  grid  in  the  spectral  domain  of  the  scattered  field 

£(m,V,  w)  ,  where: 


K,V,W  =  -— 


'W 


AA- 


N.. 


(23) 


■M 


v.v.r 


For  the  tri-linear  interpolation,  let  us  define  three  intervals  in 
the  Fourier  space:  Akx,Ak,Ak. .  Let  us  consider  a  sample 


The  major  advantage  of  this  technique  is  the  intrinsic 
possibility  of  estimate  missing  samples  when  we  choose 

A kx  .  >  A kx  .  In  this  way,  the  reconstructed  image 

shows  fewer  artifacts  and  fewer  oscillations. 

A  way  to  recover  information  from  the  missing  samples  on 
the  sparse  dataset  can  be  accomplished  by  using  the  technique 
of  Projection  on  Convex  set  (POCS).  The  basic  idea  of  POCS 
is  to  properly  weight  the  available  samples  in  a  way  that  the 
correspondent  point  spread  function  is  minimized.  The  ideal 
point  spread  function  (PSF) 

PSF(r')  =  £  (28) 

u2+v2+»2S2A0 

coincides  with  the  impulse  response  of  the  RF  tomography 
optical  system,  i.e.  has  the  shape  of  a  3D  Bessel-sinc 
function,  and  it  can  be  computed  numerically  by  using  a  3D 
FFT  algorithm  so  that  the  PSF  is  known  for  any  value  of 
r '  in  the  region  D .  For  simplicity,  let  us  consider  an 
equivalent  scaled  problem  in  which  l/,V,WG  Z  and 
AkxAkyAk_  =  1 

However,  the  actual  PSF  generated  by  the 
S  <  NXNVNZ  available  samples  is: 


point  (UyVy w)  in  the  uniform  grid  to  be  estimated,  and  a 

non-uniform  sample  point  k . =  k'  +k*  +  k~  that  has 

been  measured:  we  can  define  an  interpolation  weighting 
factor  as  follows: 

C  =  K  (C  - u ) hy  {kl,  -  v) h.  (k;}W  -  w)  (24) 
where: 


PSF(r ')  =  J  ,29) 

5=1 

where  the  weighting  factors  ps  are  identically  equal  to  unity. 
In  principle,  by  properly  weighting  the  elements  in  the 
incomplete  summation  (29),  it  may  be  possible  that  the 
weighted  actual  PSF  is  pointwise  very  similar  to  the  ideal 
PSF: 
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PSF(x\y',z')  =  PSF(x',y\z')  (30) 

In  order  to  achieve  this  result  we  invoke  the  Projection  on 
Convex  Sets  method,  which  is  based  on  successive 
approximations  of  the  actual  PSF  until  a  stopping  criterion  is 
met. 

In  the  first  iteration,  we  estimate  the  point  spread  function  by 
imposing  all  weights  to  be  ps  =  1 .  The  next  step  is  to 
compare  poiniwi.se  the  computed  PSF  with  the  ideal  PSF  (in 
the  Space  domain).  More  exactly,  we  create  a  suitable  guard 
band  region  on  the  ideal  PSF.  If  the  values  of  the  actual  PSF 
are  falling  outside  the  guard  region,  we  force  these  values  to 
be  inside  the  guard  region.  Once  the  corrections  on  the 
computed  PSF  are  done,  we  perform  a  3D  IFFT.  The  Fourier 
transform  of  the  corrected  actual  PSF  is  generally  a  function 
that  has  values  for  any  M,  V,  w  pairs.  The  next  step  is  to  set  to 
zero  all  points  in  the  Fourier  domain  that  represent  our 
missing  samples.  In  this  way,  we  are  creating  a  new  function 
defined  only  where  actual  samples  are  located. 

We  perform  a  new  computation  of  PSF  using  the  values  ps 
previously  obtained,  and  the  POCS  iteration  continues  from 
until  a  stopping  criterion  is  met,  e.g.  the  computed  PSF  lies 
completely  within  the  guard  band  region,  or  the  w-th  iteration 
of  the  computed  PSF  does  not  improve  the  approximation  of 
the  ideal  PSF  with  respect  of  the  n-lth  iteration,  i.e.  the 
process  stalls. 

When  POCS  terminates,  the  coefficient  ps  of  the  last 
iteration  are  used  for  computing  the  non-uniform  inverse 
Fourier  transform  of  the  received  electric  field,  as  follows: 


(3!) 


s=\ 


IV.  Simulations  and  Results 


We  performed  several  simulations  using  the  methods 
described  above.  The  geometry  is  depicted  in  Fig.  3,  the 
probing  frequency  is  3MHz,  and  results  are  shown  in  Figs. 
4,5,6.  Considerations  and  further  examples  will  be  shown  at 
the  time  of  the  conference. 
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Figure  3: :  Geometry  for  the  simulation  (top  view).  Transmitters  are 
represented  with  “+",  while  receivers  arc  represented  with  “X".  The  two 
black  lines  represent  the  positions  of  the  tunnels. 
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Figure  4:  reconstructed  image  using  Fourier,  no  POCS 
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Figure  5:  reconstructed  image  using  Fourier  and  POCS 


X  Barqe  [m] 

Figure  6:  reconstructed  image  using  Tikhonov 
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